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Abstract:  State-of-the-art  numerical  methods  for  solving  moving  bound¬ 
ary  problems  arising  from  multiphase  flow  and  fluid-structure  interaction 
modeling  are  reviewed.  The  emphasis  of  the  review  is  on  robust  methods 
that  do  not  require  the  mesh  to  conform  to  the  moving  boundary.  The  im¬ 
petus  for  this  review  is  the  Coastal  and  Hydraulics  Laboratory’s  mission 
in  navigation.  Accurately  predicting  the  effect  of  a  vessel  on  a  waterway 
and  the  vessel  motion  are  required  to  support  the  navigation  mission  as 
is  accurately  predicting  wave  interaction  with  coastal  and  hydraulic  struc¬ 
tures.  A  key  part  of  making  accurate  predictions  in  these  cases  is  solving 
moving  boundary  problems.  This  report  identifies  promising  directions 
for  research  and  development  of  in-house,  state-of-the-art  computational 
models. 


Disclaimer:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes.  Citation 
of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products.  All  product 
names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to  be  construed  as 
an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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ERDC/CHL  TR-09-10  iii 


Table  of  Contents 

List  of  Figures . iv 

Preface . v 

1  Introduction . 1 

2  Background . 3 

2.1  Navier-Stokes  equations  for  incompressible,  Newtonian  fluids .  3 

2.2  Moving  boundary  models .  5 

3  Current  Trends . 14 

3.1  Recent  work  on  free  boundary  modeling .  14 

3.2  Recent  work  on  multiphase  flow  modeling .  15 

3.3  Modeling  with  unstructured  tetrahedral  meshes .  16 

4  Conclusions . 17 

References . 18 


ERDC/CHL  TR-09-10  iv 


List  of  Figures 

Figure  1.  Domain  with  three  mobile  phases .  1 

Figure  2.  Parametric,  implicit,  and  characteristic  representations  of  a  surface  in  R3 .  6 

Figure  3.  Signed  distance  level  sets  for  moving  boundary  at  z  =  1/2 .  8 

Figure  4.  Level  sets  with  corners .  13 


ERDC/CHL  TR-09-10 


v 


Preface 


This  report  is  a  product  of  the  High  Fidelity  Vessel  Effects  Work  Unit 
of  the  Navigation  Systems  Research  Program  being  conducted  at  the 

U. S.  Army  Engineer  Research  and  Development  Center,  Coastal  and  Hy¬ 
draulics  Laboratory. 

The  report  was  prepared  by  Drs.  Christopher  E.  Kees  and  Matthew  W.  Far¬ 
thing  under  the  supervision  of  Mr.  Earl  V.  Edris,  Jr.,  Chief,  Hydrologic 
Systems  Branch;  by  Dr.  Tahirih  C.  Lackey  under  the  supervision  of  Dr.  Ty 

V.  Wamsley,  Chief,  Coastal  Processes  Branch;  and  by  Dr.  R.  C.  Berger 
under  the  supervision  of  Dr.  Robert  T.  McAdory,  Chief  Estuarine  En¬ 
gineering  Branch.  General  supervision  was  provided  by  Mr.  Thomas 

W.  Richardson,  Director,  CHL;  Dr.  William  D.  Martin,  Deputy  Director, 
CHL;  and  Mr.  Bruce  A.  Ebersole,  Chief  Flood  and  Storm  Protection  Divi¬ 
sion. 

Technical  advice  needed  to  complete  this  work  was  provided  by  Drs.  Stacy 
E.  Howington  and  Robert  S.  Bernard,  Coastal  and  Hydraulics  Labora¬ 
tory.  Typesetting  assistance  was  provided  by  Ms.  Sophia  K.  Vollor,  Coastal 
and  Hydraulics  Laboratory. 

Mr.  James  E.  Clausner,  Navigation  Systems  Program  Manager,  was  the 
project  manager  for  this  effort.  Mr.  W.  Jeff  Lillycrop  was  the  Technical 
Director. 

COL  Gary  E.  Johnston  was  Commander  and  Executive  Director  of  the 
Engineer  Research  and  Development  Center.  Dr.  James  E.  Houston  was 
Director. 

This  report  was  typeset  by  the  authors  with  the  HTp;X  document  prepa¬ 
ration  system.  The  report  uses  the  erdc  document  class  and  mathgifg 
fonts  package  developed  by  Dr.  Boris  Veytsman  under  the  supervision  of 
Mr.  Ryan  E.  North,  Geotechnical  and  Structures  Laboratory.  The  pack¬ 
ages  are  available  from  http://ctan.tug.org. 


ERDC/CHL  TR-09-10 


1 


1  Introduction 


Consider  the  problem  of  simulating  a  solid  (e.g.  a  vessel)  moving  in  a 
waterway  (Figure  l).  Suppose  that  the  primary  goal  of  the  simulation  is 
to  describe  the  interaction  of  the  solid’s  movement  on  the  air  and  wa¬ 
ter  in  some  sufficiently  large  box  containing  the  waterway.  The  interior 
of  the  box  is  the  domain,  and  the  boundaries  of  the  box  are  the  domain 
boundaries.  For  simplicity,  we  assume  that  the  solid  is  a  rigid  body.  It 
does  not  deform,  but  its  motion  may  depend  on  surface  forces  exerted  by 
the  fluid  and  on  gravitational  forces.  The  incompressible  Navier-Stokes 
equations  will  be  used  to  describe  both  the  water  and  air  flow,  though 
in  full-scale  applications  a  system  of  Reynolds  Averaged  Navier-Stokes 
equations  (RANS)  would  be  required.  This  problem  contains  three  mov¬ 
ing  boundaries  at  the  interfaces  between  the  three  material  types:  air- 
water,  air-solid,  and  water-solid.  The  movement  is  relative  to  the  fixed 
domain  boundaries,  and  the  interfaces  remain  sharp.  No  significant 
amounts  of  material  pass  through  the  moving  boundaries  at  the  tem¬ 
poral  and  spatial  scales  of  interest.  That  is,  the  moving  boundaries  are 
material  surfaces.  On  the  other  hand,  mechanical  and  thermal  energy 
can  and  do  pass  through  these  moving  interfaces.  The  particular  way 
that  mass  and  momentum  fluxes  are  constrained  at  the  interfaces  con¬ 
stitute  the  boundary  conditions  for  the  Navier-Stokes  equations  along 
these  moving  boundaries. 
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Domai 
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Figure  1.  Domain  with  three  mobile  phases. 


A  classic  method  for  describing  this  moving  boundary  problem  is  to  shift 
from  the  fixed  frame  of  reference  associated  with  the  box  to  a  material 
or  Lagrangian  frame  that  moves  with  the  fluids  or  solid.  Due  to  the  as- 
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sumption  that  interfaces  are  material  boundaries,  the  moving  bound¬ 
aries  in  the  fixed  frame  become  fixed  boundaries  in  the  moving  frame 
and  vice  versa.  This  approach  may  complicate  the  application  of  bound¬ 
ary  conditions  at  the  domain  boundary,  and  in  the  present  case  it  is  fur¬ 
ther  complicated  by  the  fact  that  the  moving  boundaries  are  moving  rel¬ 
ative  to  one  another.  The  transformation  from  fixed  to  moving  frame 
becomes  a  part  of  the  solution  to  the  problem-it  is  based  on  integra¬ 
tion  of  the  flow  velocity.  Instead  of  a  purely  Lagrangian  perspective, 
so-called  Arbitrary  Lagrangian-Eulerian  (ALE)  space-time  formulations 
that  blend  the  Lagrangian  perspective  near  the  moving  boundaries  with 
the  fixed  “Eulerian”  frame  away  from  moving  boundaries  are  often  pre¬ 
ferred  (e.g.  (Donea,  1982;  Hughes  and  Hulbert,  1998)).  The  ALE  formu¬ 
lation  requires  additional  decisions  about  the  blending  of  Eulerian  and 
Lagrangian  frames  that  do  not  follow  from  the  physics  alone.  These  deci¬ 
sions  form  part  of  the  solution  method. 

In  many  cases,  the  resulting  flow  dynamics  may  cause  the  material  do¬ 
mains  to  undergo  changes  in  topology.  For  instance,  breaking  waves  or 
flow  over  the  top  of  the  vessel  may  change  the  connectedness  of  the  wa¬ 
ter  and  air  domains.  The  time-dependent  space  transformation  in  the 
ALE  formulations  becomes  degenerate  during  such  topological  changes. 
Furthermore,  even  when  no  topological  changes  occur,  the  numerical 
ALE  formulation  must  be  carefully  designed  to  prevent  numerical  insta¬ 
bilities  from  polluting  the  results. 

This  review  focuses  instead  on  the  fixed  or  Eulerian  perspective  because 
of  the  difficulties  associated  with  ALE  formulations  of  the  moving  bound¬ 
ary  problem.  We  also  exclude  various  other  particle-based  numerical 
methods  that  can  be  used  to  approximate  multiphase  flow  models  (e.g. 
the  Lattice-Boltzmann  method  (Frisch  et  al.,  1986)  and  the  particle  finite 
element  method  (Idelsohn  et  al.,  2006)).  For  Eulerian  approaches,  the 
primary  difficulty  becomes  accurately  locating  the  moving  boundary  and 
applying  boundary  conditions  there.  In  the  next  section  we  provide  some 
mathematical  background  on  free  boundary  problems.  Then  we  review 
and  discuss  the  current  trends  in  numerical  methods.  We  conclude  with 
a  list  of  recommendations  for  further  numerical  modeling  work. 
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2  Background 

We  begin  by  reviewing  the  basic  mathematical  framework  of  moving 
boundary  problems  for  two-phase  flow.  First,  we  describe  the  Eulerian 
flow  model  and  then  popular  techniques  for  modeling  the  moving  bound¬ 
ary. 

2.1  Navier-Stokes  equations  for  incompressible,  Newtonian  fluids 

We  partition  the  spatial  domain,  O,  into  two  subdomains  Q,w  (water)  and 
Qa  (air).  In  each  subdomain  we  write  the  incompressible  Navier-Stokes 
equations  as: 


V  •  va  =  0  (1) 

9(P»V«)  +  v  _  (  Va  0  va  -  Oa)  -  Pag  =  0  (2) 

ot 

where  a  =  w,a  and  oa  is  the  stress  tensor,  given  by 

-  ( dVi  dvj\ 

<7,,  =  -PS,,+^—  +  — )  (3) 

We  let  r(t)  =  d£lw  f)  90a.  Note  that  both  the  subdomains  and  T(f)  are 
actually  time  dependent.  As  a  simple  example,  the  boundary  conditions 
on  r(t)  could  be 


Va-Vu,  =  0  (4) 

Pa-Pw  =  0  (5) 

The  first  condition  reflects  continuity  of  the  normal  and  tangential  com¬ 
ponents  of  the  velocity  and  the  second  continuity  of  the  pressure.  Or, 
equivalently,  these  jump  conditions  in  the  unknowns  follow  from  an  as¬ 
sumption  of  continuity  of  the  mass  and  momentum  flux.  While  it  typi¬ 
cally  makes  sense  to  assume  that  the  interface  has  no  mass  and  therefore 
that  the  mass  flux  is  continuous  at  points  on  T,  the  well-known  phenom¬ 
ena  of  surface  tension  implies  that  the  interface  can  absorb  or  impart 
momentum.  Hence,  a  more  general  set  of  jump  conditions  is 


(va  -  vj  •  n 

(oa  -  aw) ■  n 


0 

f 


(6) 

(7) 
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where  f  is  the  the  force  due  to  the  interface  and  n  is  the  outward  unit 
normal  vector  on  £lw.  Surface  tension  results  in  a  force  on  the  fluid  di¬ 
rected  entirely  in  the  normal  direction  and  typically  given  by 

f  =  Yw,aKn  (8) 

where  k  is  the  mean  curvature  of  T(f)  defined  as 

k  =  -V  •  n  (9) 

and  Yw,a  is  the  surface  tension  coefficient  for  the  air/water  interface.  For 
related  formulations  see  (Harlow  and  Welch,  1965;  Popinet  and  Zaleski, 
1999;  Tornberg  and  Enquist,  2000). 


Weak  formulation 


In  order  to  formulate  numerical  methods  and  incorporate  the  jump  con¬ 
ditions,  we  will  use  a  weak  formulation  of  the  Navier-Stokes  equation. 
For  now  we  consider  a  domain  O  and  assume  that  T  cuts  through  O  and 
divides  it  into  two  domains  Qw  and  Oa.  Multiplying  the  continuity  equa¬ 
tion  by  a  smooth  test  function  w  with  compact  support  in  Q,  integrating 
by  parts,  and  applying  the  jump  condition  on  the  normal  component  of 
the  velocity  yields 


/  vw  ■  VwdV  +  /  vww  ■  n wdS  -  I  va  ■  VwdV  +  /  vww  ■  n adS 

J  £iw  Jt  Jna  Jt 


-  /  v  •  VwdV  -  /  (va  -  vw)w  ■  n dS 
Jo.  Jt 


v • VwdV 


10 


where  we  have  defined  a  global  fluid  velocity 


v  = 


Vw  ^  ~  r 

Va  X  G 


Treating  the  momentum  equation  in  the  same  manner  yields 

'd(pv) 


10 


dt 


Pg 


wdV  -  /  (pv  0  v  -  o)  •  Vw  =  /  fwdS 
Jo  Jt 


0  (10) 
0  (11) 
0  (12) 


(13) 


(14) 
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where 


P  = 


P  = 


a  = 


Pw 

X  £  Plw 

Pa 

X  £ 

Pw 

X  £  Plw 

Pa 

X  £  Oa 

Ow 

X  £ 

Oa 

X  £  Oa 

(15) 


(16) 


(17) 


Note  that  in  general  we  only  know  that  v  •  n  is  continuous  at  points  on  T 
while  v,  p,  p,  and  p  may  be  discontinuous.  Since  the  viscosity  is  always 
positive  and  f  acts  normal  to  the  interface,  however,  v  may  be  continu¬ 
ous. 


2.2  Moving  boundary  models 

We  require  a  representation  for  the  moving  boundary,  T(f),  and  an  equa¬ 
tion  for  its  evolution.  As  noted  in  (Sethian,  2001),  there  are  three  classic 
mathematical  perspectives  on  representing  the  boundary.  First,  T(f)  can 
be  given  in  parametric  form, 

[x(£,  77,  t),  y(£,  77,  t),  z(£,  77, 0]  =  x(£,  77,  t )  (18) 

where  (£,  77)  £  R2  are  surface  coordinates.  For  example,  the  surface  de¬ 
picted  in  figure  2  can  be  parameterized  as 

x(x)  =  (f, 77,  'J'l-l2-  772)  (19) 

where  (£,  77)  £  [-1, 1]  x  [-1, 1],  Second,  T(f)  can  be  given  in  implicit  form, 

u(x,  y,  z,  t )  =  u(x,  f)  =  0  (20) 

In  the  example  above  we  could  use 

u{x,  y,z)  =  2-x2  -  y2  -  z2  (21) 

We  will  call  u  a  level  set  (LS)  function  for  the  surface  since  the  surface 
is  given  by  the  zero  LS  of  the  function  u.  Third,  instead  of  representing 
T(0,  the  subdomains  Qw  and  £2a  can  be  represented  in  set-theoretic  form 
via  their  characteristic  functions  xw( 0  and  XaiX)  defined  by 

1  for  x  £  nw 
0  for  x  £  Oa 


Xu,(x)  = 


(22) 
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Figure  2.  Parametric,  implicit,  and  characteristic  representations  of  a 

surface  in  R3. 


and  Xa  =  1  ~  Xw-  For  the  example  above,  assuming  that  water  lies  below 
air,  we  could  use 


XW  =  Xa(x)  = 


0  if  2  -  x2  -  y2  -  z2  >  0 
1  if  2  -  x2  -  y2  -  z2  <  0 


(23) 


It  is  convenient  to  consider  the  basic  numerical  methods  as  proceeding 
from  one  of  these  three  perspectives 


Level  set  methods 


First,  note  that  if  the  initial  surface  is  described  as  u(x)  =  0  then  there  is 
no  reason  why  this  u(x)  cannot  be  defined  for  all  xefl.  This  is  apparent 
in  the  3D  example  above.  Given  an  initial  condition  for  u  on  O,  we  now 
wish  to  define  equations  for  its  evolution  on  all  of  O  as  well.  Let  x(t)  be 
a  fluid  particle  path  such  that  x(0)  is  on  T(t).  It  can  be  shown  that  x(f) 
remains  on  T(f)  for  all  time.  Therefore  u(x(f),  t )  =  0.  Differentiating  with 
respect  to  t  yields 


du  ^  dx 
—  +  Vu  •  —  =  0 
dt  dt 


(24) 


Again  by  definition,  the  velocity  of  water  particles  at  the  boundary  is 
vu>  =  Furthermore  by  continuity  of  the  normal  component  of  the  ve- 
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locities,  we  know  that  on  T,  Vu  •  vw  =  Vu  •  va  so  we  can  use  either  velocity 
and  henceforth  drop  the  subscript.  Hence  u  is  related  to  v  through 


— +  Vu-v  =  0  (25) 

ot 

We  denote  this  the  non-conservative  interface  (NCI)  equation.  Note  that 
it  is  known  as  the  kinematic  boundary  condition  for  the  free  surface  T. 
We  know  from  calculus  that  n  =  Vu/||  Vu||  where  ||  •  ||  is  the  Euclidean 
norm.  The  normal  speed  is  then 


Vu 
'  || Vu 


(26) 


We  can  thus  rewrite  the  NCI  equation  in  terms  of  the  normal  speed  as  a 
Hamilton-Jacobi  interface  (HJI)  equation 


Vu||sn 


0 


(27) 


For  the  special  case  of  V  •  v  =  0  we  can  also  rewrite  the  NCI  equation  as 
a  conservative  interface  (Cl)  equation 


du 

~dt 


+  V  •  (uv)  =  0 


(28) 


Note  that  v  need  only  agree  with  the  physical  velocity  along  T,  and  in 
fact  many  implicit  representions  of  T  are  possible  that  agree  only  on  T. 
For  example,  suppose  O  is  [-1, 1]  x  [-1, 1]  and  T  is  simply  the  circle  of 
radius  1/2.  A  signed  Euclidean  distance  to  this  circle  makes  sense  for  all 
(x,  y).  The  zero  LS  of  this  signed  distance  function  is  an  implicit  repre¬ 
sentation  of  the  boundary.  In  this  simple  case  the  LS  function  is  given 
by  u(x,  y)  =  1/2  -  \/x2  +  y2.  The  level  sets  are  shown  in  Figure  3.  Note 
that  u2  =  0  is  a  different  LS  function  that  still  has  the  same  zero  LS  as 
u.  The  signed  Euclidean  distance  from  x  to  the  set  T  is  generally  a  good 
choice  for  the  LS  rebecause  it  provides  additional  geometric  information: 
its  value  gives  both  the  distance  to  the  interface  from  a  point  and  which 
side  of  the  interface  the  point  lies  on.  It  can  be  shown  that  the  signed 
distance  function  is  the  unique  vanishing  viscosity  solution  of  the  bound¬ 
ary  value  problem  for  the  eikonal  equation: 


Vu 


u 


1  for  x  e  O 
0  for  x  e  T 


(29) 

(30) 
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where  ||  ||  is  the  Euclidean  norm.  Solutions  of  the  eikonal  equation  can  in 
turn  be  characterized  as  steady  state  solutions  of  the  initial  value  prob¬ 
lem 


ut  +  sgn(u0)(||Vu||  -  1)  =  0 

u(x,0)  =  u0(x) 


(31) 

(32) 


where  Uo  is  any  LS  function  describing  T  and  sgn  is  the  signum  function 
given  by 


{-1  u  <  0 
0  u  =  0 
1  u  >  0 


(33) 


Figure  3.  Signed  distance  level  sets  for  moving 
boundary  at  z  =  1/2. 


One  additional  element  is  needed  to  make  this  description  complete:  the 
extension  of  sn  or  v  to  all  of  Q.  This  extension  is  not  unique  because  it 
need  only  match  the  normal  of  the  fluid  velocities  on  T.  Two  natural 
choices  are  the  fluid  velocities  themselves  and  a  speed  sn  that  preserves 
the  signed  distance  structure  of  u. 

The  modern  LS  method  based  on  the  HJI  equation  was  first  formulated 
in  (Osher  and  Sethian,  1988).  The  LS  method  can  lead  to  loss  of  mass 
because  no  discrete  mass  conservation  properties  can  be  derived  from 
the  fluid  and  LS  equations  alone  (even  the  Cl  formulation).  This  is  a  sig¬ 
nificant  drawback  of  the  basic  LS  method,  which  we  discuss  briefly  be¬ 
low. 
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To  complete  the  method  for  modeling  two-phase  flow  we  combine  the  LS 
representation  with  the  two-phase  Navier-Stokes  formulation.  As  an  ex¬ 
ample  we  choose  to  preserve  a  signed  distance  representation  of  the  LS 
using  the  eikonal  equation  and  use  the  NCI  equation  for  the  LS  dynam¬ 
ics.  The  complete  system  of  equations  in  weak  form  is 


la 


(|| Vu||  -  1)  wdV  =  0 


-  /  v  ■  VwdV  =  /  qpdS 

J 12  J  Tn  p 


(34) 

(35) 


la 


d(pv) 

dt 


-Pg 


where 


wdV  -  j  (pv0v-o)-Vm 

/(|  +  Vu.v|W  -  0 


/  fdS  +  /  q vdS  (36) 

It  Jtn,v 

(37) 


P  = 
P  = 

H(u)  = 


PaH(u )  +  pw  [1  -  H(u)] 
paH(u )  +  pw[l-  H{u )] 

0  u  <  0 


|  u  =  0 
1  u  >  0 


(38) 

(39) 

(40) 


and  the  boundary  and  initial  conditions  are 


u 

P 

v 

v  •  n 

(pv  ®  v  -  o)  •  n 

P 

v 

A 

u 


u  for  x  6  T 

(41) 

pd  for  x  e  r D,p 

(42) 

\d  for  x  e  Td,v 

(43) 

<?P(x,  t)  for  x  G  TN>p 

(44) 

Qd  for  x  G  TNiV 

(45) 

Po  for  x  G  O,  t  =  0 

(46) 

v0  for  x  G  O,  t  =  0 

(47) 

u0  for  x  G  O,  t  =  0 

(48) 

Note  that  equation  34  “redistances”  or  “re-initializes”  u  so  that  the  flow 
equations  are  always  able  to  depend  on  a  LS  description  that  provides 
the  signed  distance  to  the  interface.  This  system  is  typically  solved  using 
operator  splitting.  For  example,  the  solution  at  time  At  can  be  approx¬ 
imated  by  first  solving  equation  34  (redistance),  then  solving  equations 
35  and  36  (Navier-Stokes)  using  u( 0),  then  finally  equation  37  (LS)  at 
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t  =  At  using  v(A t).  A  common  variant  of  this  approach  is  to  use  smooth 
(regularized)  approximations  of  the  Heaviside  function,  H,  and  the  Dirac 
delta  function  to  evaluate  the  discontinuous  coefficients  and  surface  force 
integral: 


p  *  PaHSju )  +pw[  1  -  H£(u)] 
P  ~  paHe(u )  +pw[  1  -  H£(u)] 

fdS  ~  /  5£(u)f||Vu||dV 


in 


HM  = 


SM  = 


0 

1 

2 


0 

1 

2 

0 


u  < 

1  +  -  +  ^±2')  u  =  0 

e  n 


u  >  e 


1  cos(^) 

e  e 


u  <  -e 
u  =  0 

u  >  e 


(49) 

(50) 

(51) 


(52) 


(53) 


where  e  is  a  mesh-dependent  parameter  (e.g.  e  =  1.5Ax  is  recommended 
in  (Osher  and  Fedkiw,  2003)).  It  was  pointed  out  in  (Tornberg  and  En- 
quist,  2003)  that  these  simple  regularizations  of  the  Heaviside  and  Dirac 
functions  are  inaccurate  and  that  significantly  better  alternatives  exist. 
Furthermore,  direct  evaluation  of  the  surface  force  term  is  also  feasi¬ 
ble  as  are  related  immersed  boundary  methods  described  in  (Li  and  Ito, 
2006).  It  is  also  common  in  surface  tension  calculations  to  introduce  a 
filtering  step  to  remove  high  frequency  oscillations  from  the  surface  nor¬ 
mal  ||Vu||  or  k  (Tornberg  and  Enquist,  2000). 


We  note  that  the  primary  advantage  of  the  LS  formulation  is  that  accu¬ 
rate  and  efficient  finite  element  and  finite  volume  methods  exist  for  solv¬ 
ing  equations  34-37.  Such  methods  automatically  compute  the  viscosity 
solutions  for  the  LS  functions,  which  are  the  correct  solutions  for  cases 
in  which  the  connectivity  of  the  fluid  regions  change  (i.e.  during  break¬ 
ing  or  merging  of  the  fluid  regions). 


Before  continuing  we  consider  the  conservation  properties  of  the  basic 
LS  method.  The  mass  of  water  in  the  domain  is  given  by 


pwdV  = 


pwH(u)dV 


(54) 
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The  integral  form  of  mass  conservation  can  likewise  be  extended  to  the 
entire  domain  or  any  subdomain  as 


(55) 


Of  course,  for  incompressible  fluids  this  implies  volume  conservation  as 
well  so  that  we  could  cancel  out  pw.  Neither  a  global  nor  a  local  state¬ 
ment  of  mass  or  volume  conservation  can  be  obtained  from  any  of  the 
LS  equations  derived  above. 

Front  tracking  (parametric)  methods 

Next  we  consider  the  parametric  representation  of  the  surface,  x(£,  77,  fl¬ 
it  is  easy  to  show  from  the  implicit  representation  of  T(f)  that  ^  •  n  =  s„, 
hence  the  front  tracking  (FT)  formulation  can  be  written  as 


dx 


(56) 


where  vFT  is  any  velocity  such  that  vFT  •  n  =  sn.  Given  a  discrete  repre¬ 
sentation  of  the  x(f,  77,  f)  the  FT  equation  is  set  of  (uncoupled)  ordinary 
differential  equations.  Early  examples  of  FT  for  two-phase  flow  include 
(Glimm  et  al.,  1987;  Tryggvason  and  Unverdi,  1990;  Unverdi  and  Tryg- 
gvason,  1992).  The  primary  drawback  of  FT  is  that  for  many  realistic  ve¬ 
locities,  the  trajectories  defined  by  \FT  intersect,  leading  to  degeneracy  in 
the  parametric  representation.  These  problems  can  be  dealt  with  using 
additional,  rather  complex,  “untangling”  algorithms  based  on  additional 
physical  and  geometric  information  (Glimm  et  al.,  1988).  We  consider 
these  methods  Eulerian  methods  because  the  underlying  mesh  for  resolv¬ 
ing  physics  is  fixed,  even  though  the  Lagrangian  particle  tracking  is  a  key 
part  of  the  method.  Given  an  implementation  of  a  complete  algorithm, 
the  two-phase  model  is  similar  in  overall  design  to  the  LS  model  above, 
with  the  exception  that  the  LS  advection  step  is  replaced  by  the  FT  step 
and  the  coefficients  are  evaluated  using  the  parametric  form  of  the  sur¬ 
face. 


Volume  of  fluid  methods 


This  approach  could  be  described  as  “pseudo-mixture  theory”  because 
the  set  characteristic  functions  Xw  and  are  derived  from  the  volume 
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fraction  of  each  discrete  cell  or  element,  as  if  it  were  filled  with  a  mix¬ 
ture  of  the  two  fluids.  Instead  of  solving  the  free  boundary  problem  ex¬ 
plicitly,  the  entire  fluid  domain  is  described  by  a  single  set  of  equations 
for  the  two-fluid  mixture  with  a  particular  choice  of  “constitutive  laws” 
for  the  mixture  properties.  Let  6e(x,  t)  be  the  ratio  of  water  volume  to 
total  fluid  volume  in  a  computational  element  e.  If  initially  all  elements 
contain  either  water  or  air  then  6e  =  Xw-  We  then  define  the  density  and 
viscosity  of  the  two-fluid  mixture  as 

p  =  6pw  +  (±-6)pa  (57) 

P  =  6pw  +  (l-  6)pa  (58) 

In  a  sense  these  are  “constitutive  laws”  of  the  mixture,  though  they  are 
introduced  as  a  numerical  convenience  as  are  the  smoothed  Heaviside 
functions  in  the  LS  methods.  The  flow  equations  for  the  mixture  are  sim¬ 
ply  the  Navier-Stokes  equations  for  the  fluid  mixture,  which  necessarily 
has  variable  density  and  viscosity  (though  it  is  still  incompressible).  An 
additional  transport  equation  for  9W  represents  conservation  of  volume 
for  either  phase: 


0f  +  V  •  (0v)  =  0  (59) 

where  v  is  the  velocity  of  the  fluid  mixture.  Since  mass  conservation  is 
explicitly  enforced  one  need  only  discretize  carefully  in  order  to  conserve 
mass  in  the  fully  discrete  model.  On  the  other  hand,  Eulerian  methods 
for  the  volume  of  fluid  (VOF)  equation  tend  to  smear  the  jumps  in  6, 
which  results  in  a  loss  of  information  about  the  precise  location  of  T.  In¬ 
terface  reconstruction  and  sharpening  techniques  must  be  used  to  coun¬ 
teract  this  problem  and  to  obtain  precise  information  on  the  interface 
geometry.  Early  work  on  VOF  methods  includes  (Noh  and  Woodward, 
1976;  Hirt  and  Nichols,  1981) 


Relationships  among  interface  representations 


We  now  relate  the  FT  and  VOF  interface  representations  to  the  LS  rep¬ 
resentation.  The  simplest  relationship  is  among  solutions.  Given  a  para¬ 
metric  representation  xmp(£,  77,  t)  then  the  signed  distance  LS  function  is 
given  by  u(x,  t )  =  min^^  ||x  -  xmp(£,  77,  t)\\.  Given  a  signed  distance  LS 
function  u(x,  t)  (u  <  0  being  the  water  phase)  the  characteristic  function 
is  6  =  1  -H(u(x,  t ))  where  H  is  the  Heaviside  function  defined  in  equation 
40 
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The  methods  for  obtaining  solutions  are  also  related.  We  can  rewrite  the 
HJI  equation  as 

ut  +  Vu  •  (snn)  =  0  (60) 

This  equation  is  the  directional  derivative  along  the  curve  determined  by 
the  characteristic  equations 


Sn  n  (61) 

x(s,  0)  (62) 

which  are  equations  for  a  FT  representation.  Thus,  if  the  FT  equations 
can  be  solved  then  the  HJI  equation  can  be  as  well  via 
u(xs(t),  t )  =  u(xs( 0),  0)  and  whatever  method  is  used  to  extend  u  to  all  of 
space.  It  is  well  known  that  the  method  of  characteristics  does  not  yield 
the  solution  of  the  HJI  equations  for  large  time  because  the  character¬ 
istics  trajectories,  xs(f),  will  cross  in  many  settings  as  mentioned  above. 
On  the  other  hand,  the  theory  of  viscosity  solutions  yields  a  method  for 
extending  solutions  of  the  HJI  equation  past  the  point  when  character¬ 
istics  cross.  For  example,  far  enough  from  a  curved  boundary  the  signed 
distance  LS  function  will  form  corners,  see  Figure  4 


dxs 

dt 

Xs(0)  = 


As  we  have  seen,  the  volume  fractions  can  be  derived  directly  from  the 
LS  function  via  the  Heaviside  function.  Since  the  VOF  model  is  a  conser¬ 
vation  law,  it  can  also  be  interpreted  in  terms  of  a  method-of-characteristics 
solution  which  may  break  down  but  can  be  continued  via  the  viscosity 
solution,  similarly  to  the  HJI  equation.  The  volume  fraction  drops  some 
of  the  information  provided  in  the  LS  function,  however,  using  only  its 
sign. 
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3  Current  Trends 


Each  basic  method  or  viewpoint  has  strengths  and  weaknesses: 


•  FT  methods  give  very  accurate  representations  of  smooth  interfaces. 
They  require  complex  untangling  algorithms  to  deal  with  topological 
changes  and  are  not  strictly  mass  or  volume  conserving. 

•  LS  methods  that  deal  with  topological  changes  are  easy  to  implement. 
They  are  not  strictly  mass  or  volume  conserving  and  tend  to  be  less 
accurate  than  FT  methods. 

•  VOF  methods  are  strictly  mass  conserving.  They  can  deal  with  topo¬ 
logical  changes,  but  they  tend  to  be  the  least  accurate  of  the  three 
methods. 


The  following  advice  from  (Sethian,  2001) 


Good  numerics  is  ultimately  about  getting  things  to  work;  the 
slavish  and  blind  devotion  to  one  approach  above  all  others  is 
usually  a  sign  of  unfamiliarity  with  the  range  of  troubles  and  chal¬ 
lenges  presented  by  real  applications. 


Much  of  the  recent  work  on  Eulerian  methods  reflects  attempts  at  com¬ 
bining  aspects  of  all  three  methods  in  order  to  obtain  the  best  properties 
of  all  the  methods. 

3.1  Recent  work  on  free  boundary  modeling 

One  approach  to  improve  the  mass  conservation  properties  of  LS  meth¬ 
ods  is  the  conservative  LS/VOF  (CLSVOF)  method  (Sussman  and  Puck¬ 
ett,  2000).  In  this  approach  both  the  VOF  and  LS  representations  are 
computed  and  the  VOF  function  is  used  to  obtain  the  correct  mass.  The 
CLSVOF  method  is  currently  being  modified  and  extended  to  problems 
with  surface  tension  (Sussman,  2003;  Sussman  et  al.,  2007).  A  related 
approach  simply  uses  a  LS  function  that  mimics  a  smoothed  Heaviside 
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function  rather  than  a  signed  distance  function  in  order  to  mimic  the 
shape  of  the  set  characteristic  function  (Olsson  and  Kreiss,  2005). 

A  different  approach  to  maintaining  mass  conservation  and  improving 
the  resolution  of  LS  methods  is  the  particle-LS  method  presented  in  (En¬ 
right  et  al.,  2002).  This  approach  seeds  both  sides  of  the  interface  with 
“fluid”  particles  that  are  then  used  to  correct  the  interface.  This  method 
improves  the  mass  conservation  properties  of  the  method.  Lagrangian 
tracking  ideas  are  also  being  explored  to  improve  the  accuracy  of  low  or¬ 
der  methods  for  the  LS  advection  (Enright  et  al.,  2005). 

Much  of  the  research  on  VOF  methods  has  focused  on  obtaining  higher 
accuracy  and  better  representations  of  the  interface  geometry.  The  orig¬ 
inal  piecewise  linear  interface  reconstruction  technique  (PLIC)  (Youngs, 
1982)  has  been  improved  upon  using  parabolic  (PROST)  (Renardy  and 
Renardy,  2002)  and  least  squares  (Puckett  et  al.,  1997)  techniques.  A 
nice  set  of  comparisons  presented  in  (Gerlach  et  al.,  2006)  indicated  that 
the  CLSVOF  was  generally  superior  than  most  state-of-the-art  recon¬ 
struction  techniques  based  purely  on  VOF.  Other  notable  comparisons 
of  reconstruction  algorithms  include  (Rider  and  Kothe,  1995;  Pilliod  and 
Pucket,  1997;  Rider  and  Kothe,  1998;  Scardovelli  and  Zaleski,  2003) 

The  FT  research  has  also  continued  to  mature  so  that  three-dimensional 
simulations  of  problems  with  complex  topological  changes  are  now  pos¬ 
sible.  A  recent  set  of  comparisons  and  a  freely  available  library  of  FT 
tools  are  described  in  (Du  et  al.,  2006). 

Other  recent  directions  in  LS  research  includes  the  use  of  various  meth¬ 
ods  for  enriching  the  finite  element  space  locally  (Coppola-Owen  and  Co- 
dina,  2005;  Chessa  and  Belytschko,  2003).  The  particle  finite  element 
methods  of  (Idelsohn  et  al.,  2006)  is  a  promising  new  Lagrangian  tech¬ 
nique  that  may  challenge  Eulerian  methods  even  with  respect  to  robust¬ 
ness. 

3.2  Recent  work  on  multiphase  flow  modeling 

Recently,  (Kadioglu  et  al.,  2005)  used  LS  methods  for  modeling  under¬ 
water  explosions  requiring  both  compressible  and  incompressible  flow 
over  a  wide  range  of  Mach  number.  Of  particular  interest  to  inland  nav¬ 
igation  is  (Price  and  Chen,  2006),  which  considers  a  standard  LS  for¬ 
mulation  of  two-phase  flow  for  free  surface  waves.  Reynolds-averaged 
Navier-Stokes  equations  are  used  to  model  the  fluid  flow  while  a  smoothed 
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Heaviside  function  is  used  for  material  properties  (surface  tension  is  ne¬ 
glected).  The  signed  distance  property  is  preserved  by  time  marching  the 
eikonal  equation  to  equilibrium.  Problems  considered  include  oscillating 
flow  in  a  tank  (sloshing),  a  dam  break  problem,  and  a  ship  hull  moving 
in  calm  water.  Another  interesting  application  is  the  physics-based  sim¬ 
ulations  described  in  (Irving  et  al.,  2006),  which  combines  some  aspects 
of  shallow-water  modeling  with  three-phase  Navier-Stokes/solid  body 
modeling. 

3.3  Modeling  with  unstructured  tetrahedral  meshes 

Most  of  the  the  methods  above  are  compatible  with  unstructured  meshes, 
but  much  of  the  work  has  been  carried  out  on  Cartesian  meshes.  The  LS 
method  was  extended  to  triangular  and  tetrahedral  meshes  in  (Barth  and 
Sethian,  1998).  That  method  was  extended  to  air/water  flow  in  porous 
media  in  (Holm  and  Langtangen,  1999).  Recently  (Nagrath  et  al.,  2005) 
use  stabilized  finite  element  for  two-phase  flow  using  LSs  with  a  smoothed 
Heaviside  function  for  material  properties.  A  discontinuous  Galerkin 
method  appropriate  for  unstructured  meshes  is  given  in  (Marchandise 
et  al.,  2006).  VOF  methods  have  also  been  extended  to  triangular  and 
tetrahedral  meshes  (Aliabadi  et  al.,  2003;  Tezduyar  et  al.,  1998). 
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4  Conclusions 


We  have  reviewed  a  range  of  methods  for  accurately  and  efficiently  solv¬ 
ing  two-  and  three-phase  boundary  problems.  There  is  as  yet  no  consen¬ 
sus  on  a  best  method.  Instead,  there  is  still  much  activity  in  the  area  of 
computational  physics  and  engineering  directed  at  deriving  improved  al¬ 
gorithms  by  combining  mass/volume  conservation,  particle  tracking,  and 
LS  methods.  Given  that  state  of  affairs  we  conclude  that 


1.  A  mass  conservative  hybrid  method  organized  around  a  core  LS  im¬ 
plementation  is  the  most  practical  choice  for  initial  development,  hav¬ 
ing  a  good  mix  of  robustness,  accuracy,  and  efficiency.  Additionally, 
the  tools  will  have  a  wide  range  of  applicability  outside  of  two-  and 
three-phase  flow. 

2.  Including  Lagrangian  (particle  path)  information  has  consistently 
been  a  route  toward  higher  accuracy  results  so  that  some  form  of  par¬ 
ticle  tracking  will  likely  also  be  involved  in  an  effective  algorithm.  We 
should  begin  experimenting  with  both  hybrid  particle-LS  methods 
and  the  more  complex  open  source  FT  implementations. 

3.  Completion  of  a  two-phase  test  set  (e.g.  (Rider  and  Kothe,  1995;  Du 
et  al.,  2006))  should  be  the  highest  priority  for  the  next  phase  of  the 
project.  A  test  set  will  allow  quantitative  studies  of  accuracy  and  ef¬ 
ficiency  on  problems  of  concern  to  inland  navigation  as  well  as  guide 
laboratory  and  field  data  collection  for  model  validation. 
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